1 Objetivo:

Realizar a modelagem geoestatística da profundidade do solum e da altura total das árvores de um povoamento florestal de Pinus taeda

2 Metodologia

2.1 Caracterização da área de estudo

Os dados que utilizei para este trabalho são provenientes de um povoamento florestal de 108 ha de Pinus taeda L. de 29 anos. A área pertence ao município de Campo Belo do Sul, região serrana do Estado de Santa Catarina, Brasil. O clima é do tipo Cfb, mesotérmico, subtropical úmido e com precipitação média de 1.647 mm com chuvas bem distribuídas no ano. A geologia regional é constituída por uma sequência vulcânica de rochas ácidas da Formação Serra Geral, com predomínio de riodacito.

Na área foram identificados Neossolos Litólicos e Neossolos Regolíticos, Cambissolos Háplicos e Cambissolos Húmicos, Latossolos Vermelhos e Gleissolos Melânicos, os quais representam a topossequência clássica da área de estudo (Figura 1c). Conforme as tendências observadas no campo, em locais de relevo plano ou suavemente ondulado com boa drenagem estão solos profundos com sequência de horizontes A-Bw (Latossolos), em condições de má drenagem, ocorrem solos com a sequência de horizonte A-Cg (Gleissolos). Nos relevos ondulado ou fortemente ondulado, predominam solos rasos com a sequência de horizontes A ou A-Bi (Neossolos e Cambissolos).

1 Localização da área de estudo no município de Campo Belo do Sul, Estado de Santa Catarina (SC), Brasil (a) e área ampliada com modelo digital de elevação e distribuição dos pontos amostrais (b). Relação da classe de solo e altura das árvores de uma Topossequência típica da área de estudo (c).

Figure 2.1: 1 Localização da área de estudo no município de Campo Belo do Sul, Estado de Santa Catarina (SC), Brasil (a) e área ampliada com modelo digital de elevação e distribuição dos pontos amostrais (b). Relação da classe de solo e altura das árvores de uma Topossequência típica da área de estudo (c).

2.2 Amostragem

A coleta de dados foi realizada em 102 pontos (Figura 1 b) alocados pelo algoritmo Hipercubo Latino Condicionado (Conditioned Latin Hipercube Sampling – cLHS) através da função clhs::clhs(). Foram consideradas como variáveis condicionantes da amostragem a elevação, profundidade do vale, índice de umidade tipográfico, nivel base da rede de drenagem e declividade as quais, juntas, explicaram aproximadamente 86% da variância topográfica, identificada através da Análise de Componentes Principais (ACP).

Em cada ponto amostral foi medimos a profundidade do solum (PF) e a altura total das árvores (h). Consideramos solum a espessura máxima do solo onde as raízes podem se desenvolver sem impedimentos físicos para penetração livre das raízes. Os fatores limitantes considerados foram o lençol freático elevado e o contato com rocha consolidada (contato lítico) com ou sem fissuras.

Para representar a altura das árvores utilizamos o valor médio dos quatro indivíduos de Pinus taeda mais próximos de cada ponto amostral. Portanto, considerando o espaçamento existente entre os indíviuos, cada ponto de amostragem representa uma área de aproximadamente 100 m² no campo.

Os dados foram armazenados no objeto pontos definido como SpatialPolygonsDataFrame com as coordenadas projetadas em WGS84.

pontos <- read.csv('../data/GateadosDados.csv', dec = ".", sep= ";", stringsAsFactors = FALSE)
sp :: coordinates(pontos) <- c('X' , 'Y')
wgs84utm22s <- sp::CRS('+proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
sp :: proj4string(pontos) <- wgs84utm22s
pontos <- sp :: spTransform(pontos, wgs84utm22s)

2.3 Modelo geoestatístico

Para isso foi utilizado o modelo linear misto de variação espacial \(y(\boldsymbol{s}_i)\), denotado por \[Y(\boldsymbol{s}_i) = Z(\boldsymbol{s}_i) + \varepsilon(\boldsymbol{s}_i) = \boldsymbol{x}(\boldsymbol{s}_i)^\text{T}\boldsymbol{\beta} + B(\boldsymbol{s}_i) + \varepsilon(\boldsymbol{s}_i)\] Segundo esse modelo, os dados geoestatísticos são uma realização de um campo aleatório \(Y(\boldsymbol{s}_i)\) que podem ser descritos como a combinação aditiva de efeitos fixos, efeitos aleatórios e erro aleatório independente.

\(Z(\boldsymbol{s}_i)\) é o chamado sinal e possui dois componentes. O primeiro, \(\boldsymbol{x}(\boldsymbol{s}_i)^\text{T}\boldsymbol{\beta}\) representa os efeitos fixos, ou seja, a tendência espacial de origem desterminística, especificamente, a relação de dependência (causa e efeito) entre as covariáveis ambientais e a variável ambiental.

O segundo componente do sinal, \(B(\boldsymbol{s}_i)\), representa os efeitos aleatórios, especificamente, um campo aleatório Gaussiano estacionário não-observável, também chamado de variável latente. Esse campo aleatório Gaussiano é inteiramente descrito por sua função de média e função de covariância.

Por fim, \(\varepsilon(\boldsymbol{s}_i)\) é um erro (ou ruído) independente e identicamente distribuído (iid), descrito por uma distribuição Gaussina de probabilidade, cujo parâmetro desconhecido de escala é \(\tau\). Tradicionalmente, \(\tau^2\) é chamado de efeito pepita ou ainda de variância residual.

2.3.1 Efeito fixo

Defini o \(\boldsymbol{x}(\boldsymbol{s}_i)^\text{T}\boldsymbol{\beta}\) como a profundidade do solum predita via floresta aleatória.

As covariáveis utilizadas na predição via floresta aleatória foram os atributos do terreno elevação (ELEV), declividade (DECLI), índice de umidade topográfico (TWI) e profundidade do vale (VD) e a informação do índice de produtividade ou índice de Sítio (Sitio).

2.3.1.1 Covariáveis topográficas

As covariáveis topográfricas foram derivadas do modelo digital de elevação disponibilizado pelo Governo do Estado de SC - Secretaria de Estado do Desenvolvimento Econômico Sustentável, proveniente do Levantamento Aerofotogramétrico em 2010. Os dados, disponibilizados com resolução de 1 metro, foram reamostrados para resolução espacial de 10 metros utilizando a ferramenta reamostragem no software SAGA GIS.

Os planos de informação foram importados pela ferramenta raster::rastere armazenados em objetos no formato RasterLayer.

DECLI <- raster::raster("../data/Covars/DECLI.tif")
ELEV <- raster::raster("../data/Covars/ELEV.tif")
TWI <- raster::raster("../data/Covars/TWI.tif")
VD <- raster::raster("../data/Covars/VD.tif")

Figure 2.2: 2 Covariáveis topográficas utilizadas como preditoras no modelo floresta aleatória

A partir da função raster::extract foram extraídos os valores de cada objeto raster nas localizações de cada observação contidas no objeto espacial pontos.

pontos$DECLI <- raster::extract(DECLI, pontos)
pontos$ELEV <- raster::extract(ELEV, pontos)
pontos$TWI <- raster::extract(TWI, pontos)
pontos$VD <- raster::extract(VD, pontos)

2.3.1.2 Covariável de produção

Além dos atributos de terreno, foi utilizado um índice de produtividade disponibilizado pela empresa como variável preditora do modelo floresta aleatória.

Esses índices são utilizados para ordenamento da produção. A área possui 12 parcelas fixas de inventário contínuo (PIC) de 500 metros quadrados cada que são utilizados na estimativa da produtividade local. Cada PIC é classificada em função da média de altura das 100 árvores de maior perÍmetro basal da parcela, expressa como altura dominante (Hdom). Em função da Hdom e seus incrimentos anuais é estabelecido o índices de sítio (Sitio). Esse índice é então atribuido à todo o talhão.

As parcelas da área foram classificação em 4 níveis, em que 1 corresponde ao sítio com melhor produtividade. Esses níveis de produtividade foram atribuidos aos polígonos e, naquelas em que não existe parcela de PIC, o valor de sítio foi estimado com base na altura das árvores amostradas á campo.

Utilizei a função raster::shapefile para carregar o polígono com informações das sítio - armazenado no objeto ProdutividadePistola. A função sp::spTransform foi usada para projetar as coordenadas original no plano cartesiano (UTM) e a função raster::extract para extraír os valores de cada objeto raster nas localizações de cada observação contidas no objeto espacial pontos.

Sitio <- raster::raster("../data/Covars/Sitio.tif")
sp::proj4string(Sitio) <- wgs84utm22s
pontos$Sitio <- raster::extract(Sitio, pontos) %>%  as.factor()

ProdutividadePistola <- 
 raster::shapefile('../data/Produtividade/ProdutividadePistola.shp') %>%
 sp::spTransform(wgs84utm22s)
ProdutividadePistola$Sitio <- as.factor(ProdutividadePistola$Sitio)

sp::spplot(
  ProdutividadePistola, scales = list(draw = T),
  main = 'Classificação dos sítios e  localização dos pontos') +
  lattice::xyplot(Y ~ X, data = as.data.frame(pontos@coords),
                  pch = 20, col = 'red', lwd = 2, cex = 1.5) %>%
latticeExtra::as.layer()

Figure 2.3: 3 Classificação dos sítios e localização dos pontos na área de estudo

2.3.1.3 Predição via floresta aleatória

foi construido um modelo de predição utilizado para através da função carret::train pelo método floresta aleatório method = "rf", utilizando como covariávies preditoras as variáveis topográficas e de produtividade.

rf_fit <- caret::train((PFd ~ DECLI + ELEV + TWI + VD + Sitio), data = pontos@data, method = "rf", tuneLength = 1, ntree = 1000, importance = TRUE, na.action = na.omit, trControl = trainControl("LOOCV"))
rf_fit

## Random Forest 
## 
## 94 samples
##  5 predictor
## 
## No pre-processing
## Resampling: Leave-One-Out Cross-Validation 
## Summary of sample sizes: 93, 93, 93, 93, 93, 93, ... 
## Resampling results:
## 
##   RMSE      Rsquared  
##   3.588665  0.03069723
## 
## Tuning parameter 'mtry' was held constant at a value of 2

A Figura @ref(fig:lm-residuos) mostra o ajustade dos dados de profundidade do solum preditos e observaados.

par <- par(mfrow = c(2,2))
pontos@data$rf <- rf_fit$finalModel$predicted
lm(PFd ~ rf, data = pontos) %>% plot()

Figure 2.4: 4 Ajuste de dados preditos e observados

Predição espacial

beginCluster()

## 4 cores detected, using 3

prediction <- 
  clusterR(brick(DECLI, ELEV, TWI, VD, Sitio), raster::predict, 
           args = list(model = rf_fit, type = "raw", index = 1))
endCluster()
plot(prediction)

pontos$PFprediction <- raster::extract(prediction,  pontos)

3 Para utilizar esse modelo, supus que os dados são uma realização de um campo aleatório com distribuição normal que podem ser descritos como a combinação aditiva de efeitos fixos, efeitos aleatórios e erro aleatório independente.

3.0.1 Ajuste do variograma

distmax <-dist(pontos@coords) %>% max() /3

limites <- seq(0, distmax, length.out = 20)
vario<- georob::sample.variogram(PFd ~ PFprediction,
    data= pontos, locations = ~ X + Y, lag.dist.def = limites, estimator = "matheron") %>%
  plot(ylab = 'Seminvariância', xlab = 'Distância de separação (m)', annotate.npairs = TRUE)

lags <- seq(0, distmax, length.out = 12)
georob::sample.variogram(
 PFd ~ PFprediction, data = pontos, locations = ~ X + Y, lag.dist.def = lags,
  xy.angle.def = c(0, 22.5, 67.5, 112.5, 157.5, 180)) %>% 
    plot(type = "b", ylab = 'Semivariância', xlab = 'Distância de separação (km)')

Avaliada a evolução da semivariância nas diferentes direções, há evidência da existência de estruturas de autocorrelação espacial dependentes da direção, então o processo espacial é anisotrópico. Porém, assumi a isotropia do processo espacial e o semivariograma amostral gerado foi independente da direção.

Ajustei ao variograma amostral um modelo exponencial do variograma usando o método dos quadrados mínimos não-lineares ponderados, com ponderação definida conforme o método de Cressie. O processo de estimativa dos parâmetros do modelo exponencial do variograma foi conduzido via otimização usando a função stats::optim(method = “BFGS”).

produzidos pelo otimizador a cada 10 iterações:

vario_fit <- 
  georob::fit.variogram.model(
  vario, variogram.model = 'RMexp', param = c(variance = 10, nugget = 5, scale = 70), weighting.method = "cressie", method = "BFGS")


summary(vario_fit)

## 
## Call:georob::fit.variogram.model(sv = vario, variogram.model = "RMexp", 
##     param = c(variance = 10, nugget = 5, scale = 70), weighting.method = "cressie", 
##     method = "BFGS")
## 
## Convergence in 72 function and 29 Jacobian/gradient evaluations
## 
## Residual Sum of Squares: 63.64244 
## 
## Residuals (epsilon):
##     Min      1Q  Median      3Q     Max 
## -4.3131 -1.8347 -0.5159  1.1462  3.6763 
## 
## Variogram:  RMexp 
##                Estimate    Lower Upper
## variance       12.69262  9.10587 17.69
## snugget(fixed)  0.00000       NA    NA
## nugget          1.33044  0.05434 32.58
## scale          70.37501 50.43150 98.20

plot(vario, xlab = 'Distância de separação (m)', ylab = 'Semivariância')
lines(vario_fit, col = "red", lty = 'dashed')

Na geoestatística clássica, a variância não explicada é expressa em um único parâmetro nugget. Porém, no modelo linear misto podemos separar o parâmetro nugget em dois componentes: a variãncia devida aos erros de medida, modelada pelo parâmetro nugget e variância devida à variação espacial em pequena escala, modelada pelo parâmetro snugget. Considerando que durante a obteção dos dados de campo as informações de PF os valores foram arredondadas pelas em decímetros, considerei que a variância do erro de medida é igual a 0.5, ou seja, 62% do parâmetro nugget. A variância restante foi atribuída a variação espacial não auto-correlacionada espacialmente - não capturada pelo plano amostral snugget.

vario_fit_error <- georob::georob(
   PFd ~ PFprediction, pontos, locations = ~ X + Y, variogram.model = 'RMexp', 
 param = c(variance = vario_fit$variogram.object[[1]]$param[['variance']], 
           nugget = vario_fit$variogram.object[[1]]$param[['nugget']]*0.66,
           snugget = vario_fit$variogram.object[[1]]$param[['nugget']]*0.33,
           scale = vario_fit$variogram.object[[1]]$param[['scale']]),
 fit.param = georob::default.fit.param(nugget = FALSE, snugget = FALSE),
 tuning.psi = 1000, control = georob::control.georob(initial.fixef = 'lm'))
summary(vario_fit_error)

## 
## Call:georob::georob(formula = PFd ~ PFprediction, data = pontos, locations = ~X + 
##     Y, variogram.model = "RMexp", param = c(variance = vario_fit$variogram.object[[1]]$param[["variance"]], 
##     nugget = vario_fit$variogram.object[[1]]$param[["nugget"]] * 
##         0.66, snugget = vario_fit$variogram.object[[1]]$param[["nugget"]] * 
##         0.33, scale = vario_fit$variogram.object[[1]]$param[["scale"]]), 
##     fit.param = georob::default.fit.param(nugget = FALSE, snugget = FALSE), 
##     tuning.psi = 1000, control = georob::control.georob(initial.fixef = "lm"))
## 
## Tuning constant:  1000 
## 
## Convergence in 8 function and 7 Jacobian/gradient evaluations
## 
## Estimating equations (gradient)
## 
##                                   xi         scale
##   Gradient           :  1.392651e-03  3.466821e-03
## 
## Maximized restricted log-likelihood: -161.8979 
## 
## Predicted latent variable (B):
##      Min       1Q   Median       3Q      Max 
## -1.88319 -0.44060 -0.03349  0.45975  1.44721 
## 
## Residuals (epsilon):
##      Min       1Q   Median       3Q      Max 
## -1.82465 -0.44576 -0.03924  0.45241  1.41115 
## 
## Standardized residuals:
##      Min       1Q   Median       3Q      Max 
## -2.78661 -0.68668 -0.05961  0.69378  2.15695 
## 
## 
## Gaussian REML estimates
## 
## Variogram:  RMexp 
##                Estimate   Lower   Upper
## variance         0.4621  0.1527   1.398
## snugget(fixed)   0.4390      NA      NA
## nugget(fixed)    0.8781      NA      NA
## scale           22.2009  2.4363 202.304
## 
## 
## Fixed effects coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -4.00949    0.43847  -9.144 1.44e-14 ***
## PFprediction  1.67512    0.07029  23.830  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error (sqrt(nugget)): 0.9371 
## 
## Robustness weights: 
##  All 94 weights are ~= 1.

plot(vario)
lines(vario_fit_error)

CRIAR O GRID: Suporte de predição pontual

grid <- sp::spsample(ProdutividadePistola, 10000, type = 'regular')
colnames(grid@coords) <- colnames(pontos@coords)

grid$PFprediction <- raster::extract(prediction, grid)
grid

## class       : SpatialPointsDataFrame 
## features    : 9988 
## extent      : 516809.1, 518847.7, 6908630, 6909819  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0 
## variables   : 1
## names       :     PFprediction 
## min values  :  2.2076320473944 
## max values  : 9.17872378549943

grid <- 
  sp::SpatialPointsDataFrame(
    coords = grid@coords, 
    data = data.frame(grid),
    proj4string = grid@proj4string)
colnames(grid@coords) <- colnames(pontos@coords)

pred_ponto <- predict(vario_fit_error, newdata = grid, control = georob::control.predict.georob(extended.output = TRUE))
sp::gridded(pred_ponto) <- TRUE

## Warning in points2grid(points, tolerance, round): grid has empty column/
## rows in dimension 1

sp::spplot(pred_ponto)

Argumento type: fazer a predição do sinal “signal” z, “response” se for y ou efeitos fixos “trend” se quiser predizer erro de medida, y = z + e, então response - signal = erro se tendência - predito = componente aleatório

Usar control.predict.georob(extended.output = TRUE)

pred_ponto <- predict(
  vario_fit_error, newdata = grid, type= "response", signif = 0.95,
  control = georob::control.predict.georob(extended.output = TRUE))
sp::gridded(pred_ponto) <- TRUE

## Warning in points2grid(points, tolerance, round): grid has empty column/
## rows in dimension 1

str(pred_ponto)

## Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
##   ..@ data       :'data.frame':  9988 obs. of  8 variables:
##   .. ..$ pred           : num [1:9988] 4.29 1.13 4.57 4.48 4.47 ...
##   .. ..$ se             : num [1:9988] 1.34 1.36 1.34 1.34 1.34 ...
##   .. ..$ lower          : num [1:9988] 1.66 -1.53 1.94 1.84 1.84 ...
##   .. ..$ upper          : num [1:9988] 6.92 3.78 7.2 7.11 7.1 ...
##   .. ..$ trend          : num [1:9988] 4.29 1.14 4.57 4.47 4.47 ...
##   .. ..$ var.pred       : num [1:9988] 0.0243 0.0599 0.0229 0.0234 0.0234 ...
##   .. ..$ cov.pred.target: num [1:9988] 0.000108 0.001447 0.000154 0.00017 0.000189 ...
##   .. ..$ var.target     : num [1:9988] 1.78 1.78 1.78 1.78 1.78 ...
##   .. ..- attr(*, "variogram.object")=List of 1
##   .. .. ..$ :List of 9
##   .. .. .. ..$ variogram.model: chr "RMexp"
##   .. .. .. ..$ param          : Named num [1:4] 0.462 0.439 0.878 22.201
##   .. .. .. .. ..- attr(*, "names")= chr [1:4] "variance" "snugget" "nugget" "scale"
##   .. .. .. ..$ fit.param      : Named logi [1:4] TRUE FALSE FALSE TRUE
##   .. .. .. .. ..- attr(*, "names")= chr [1:4] "variance" "snugget" "nugget" "scale"
##   .. .. .. ..$ isotropic      : logi TRUE
##   .. .. .. ..$ aniso          : Named num [1:5] 1 1 90 90 0
##   .. .. .. .. ..- attr(*, "names")= chr [1:5] "f1" "f2" "omega" "phi" ...
##   .. .. .. ..$ fit.aniso      : Named logi [1:5] FALSE FALSE FALSE FALSE FALSE
##   .. .. .. .. ..- attr(*, "names")= chr [1:5] "f1" "f2" "omega" "phi" ...
##   .. .. .. ..$ sincos         :List of 6
##   .. .. .. .. ..$ co: num 6.12e-17
##   .. .. .. .. ..$ so: num 1
##   .. .. .. .. ..$ cp: num 6.12e-17
##   .. .. .. .. ..$ sp: num 1
##   .. .. .. .. ..$ cz: num 1
##   .. .. .. .. ..$ sz: num 0
##   .. .. .. ..$ rotmat         : num [1:2, 1:2] 1.00 -6.12e-17 6.12e-17 1.00
##   .. .. .. ..$ sclmat         : Named num [1:2] 1 1
##   .. .. .. .. ..- attr(*, "names")= chr [1:2] "" "f1"
##   .. ..- attr(*, "psi.func")= chr "logistic"
##   .. ..- attr(*, "tuning.psi")= num 1000
##   .. ..- attr(*, "type")= chr "response"
##   ..@ coords.nrs : num(0) 
##   ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
##   .. .. ..@ cellcentre.offset: Named num [1:2] 516809 6908630
##   .. .. .. ..- attr(*, "names")= chr [1:2] "X" "Y"
##   .. .. ..@ cellsize         : Named num [1:2] 8.68 8.68
##   .. .. .. ..- attr(*, "names")= chr [1:2] "X" "Y"
##   .. .. ..@ cells.dim        : Named int [1:2] 236 138
##   .. .. .. ..- attr(*, "names")= chr [1:2] "X" "Y"
##   ..@ grid.index : int [1:9988] 32417 32140 32180 32181 32182 31904 31905 31906 31944 31945 ...
##   ..@ coords     : num [1:9988, 1:2] 517538 517182 517529 517538 517546 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:9988] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "X" "Y"
##   ..@ bbox       : num [1:2, 1:2] 516805 6908626 518852 6909823
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "X" "Y"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

spplot(pred_ponto)

at <- pred_ponto@data[, c("pred", "lower", "upper")] %>% range()
sp::spplot(pred_ponto, zcol = c("lower", "pred", "upper"), main = "prediction")

sp::spplot(pred_ponto, zcol = 'se', main = "SE")

validacao <- georob::cv(vario_fit_error, nset = 93)
summary(validacao)

## 
## Statistics of cross-validation prediction errors
##        me       mede       rmse       made        qne       msse  
## -0.004617  -0.079438   1.347830   1.394152   1.364190   1.024841  
##    medsse       crps  
##  0.486049   0.766357

1 - sum((validacao$pred$data - validacao$pred$pred)^2) / sum((validacao$pred$data - mean(validacao$pred$data))^2)

## [1] 0.8569975

plot(validacao)

Parte II - Altura das árvores

n = 96 observações para calibração + predição + validação (108 ha) Validação: validação cruzada (leave-one-out)

Modelo linear misto de variação espacial (georob)

Para utilizar esse modelo, supus que os dados são uma realização de um campo aleatório com distribuição normal que podem ser descritos como a combinação aditiva de efeitos fixos, efeitos aleatórios e erro aleatório independente. Assim, considerei como efeitos fixos da variação da h a predição da PF via floresta aleatória, atributos de terreno oriundos do MDE e os índices de sítio.

Ajuste do variograma:

distmax <-dist(pontos@coords) %>% max() /2

limites <- seq(0, distmax, length.out = 20)
vario<- georob::sample.variogram(h ~ DECLI,
    data= pontos, locations = ~ X + Y, lag.dist.def = limites, estimator = "matheron") %>%
  plot(ylab = 'Seminvariância', xlab = 'Distância de separação (m)', annotate.npairs = TRUE)

lags <- seq(0, distmax, length.out = 12)
georob::sample.variogram(
 h ~ DECLI, data = pontos, locations = ~ X + Y, lag.dist.def = lags,
  xy.angle.def = c(0, 22.5, 67.5, 112.5, 157.5, 180)) %>% 
    plot(type = "b", ylab = 'Semivariância', xlab = 'Distância de separação (km)')